Nonlinear diffusion model for Rayleigh- Taylor mixing 
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The complex evolution of turbulent mixing in Rayleigh- Taylor convection is studied in terms 
of eddy diffusiviy models for the mean temperature profile. It is found that a non-linear model, 
derived within the general framework of Prandtl mixing theory, reproduces accurately the evolution 
of turbulent profiles obtained from numerical simulations. Our model allows to give very precise 
predictions for the turbulent heat flux and for the Nusselt number in the ultimate state regime of 
thermal convection. 
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Turbulent thermal convection is one of the most im- 
portant manifestations of turbulence. It appears in many 
natural phenomena, from heat transport in stars to at- 
mosphere and oceanic mixing, and it also plays a funda- 
mental role in many technological applications [![. 

This Letter is devoted to the study of turbulent convec- 
tion in the Rayleigh- Taylor (RT) setup, a paradigmatic 
configuration in which a heavy layer of fluid is placed on 
the top of a light layer. Gravitational instability at the in- 
terface of the two layers leads to a turbulent mixing zone 
which grows in time at the expenses of available potential 
energy . Specific applications of RT convection range 
from cloud formation ]^|, to supernova explosion 0, Q 
and solar corona heating [f|. Because of the absence of 
boundaries, the phenomenology of RT turbulence results 
simpler than other convective systems where the thermal 
forcing is provided by walls, such as the Rayleigh-Benard 
configuration. 

Recent theoretical work confirmed by numerical 
simulations (EMI, predicts for RT turbulence at small 
scales a turbulent cascade with Kolmogorov-Obukhov 
scaling (Bolgiano scaling in two dimensions). Here we 
concentrate on large scale features of RT turbulence. We 
propose a simple closure scheme based on the general 
framework of Prandtl mixing length theory and leading 
to a nonlinear diffusion model for temperature concen- 
tration. Our closure reproduces with high accuracy the 
spatial-temporal evolution of the mean temperature pro- 
file and allows to derive a prediction for the scaling law 
of Nu versus Ra which fits perfectly data obtained from 
direct numerical simulations. 

The equation of motion for the incompressible veloc- 
ity field v (V ■ v = 0) and temperature field T in the 
Boussinesq approximation is 



v • Vv = 

- V • VT : 
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where f3 is the thermal expansion coefficient, v the 
kinematic viscosity, n the thermal diffusivity and g = 
(0, 0, —g) is the gravitational acceleration. 

The initial condition (at t = 0) is a layer of cooler 
(heavier) fluid on the top of a hotter (lighter) layer at 



rest, i.e. v(x, 0) = and T(x, 0) = — (#o/2)sgn(z) where 
9q is the initial temperature jump which fixes the Atwood 
number A = (l/2)/3#o (T = is the reference mean tem- 
perature). This configuration is unstable and after the 
linear instability phase, the system develops a turbulent 
mixing zone which grows in time starting from the plane 
z = 0. An example of the turbulent temperature field ob- 
tained from high resolution direct numerical simulations 
of (pSD is shown in Fig. [TJ 




FIG. 1: (color online). Snapshot of a (x,z) section of the 
temperature field for Rayleigh- Taylor turbulence numerical 
simulation. White (black) represents hot, light (cold, heavy) 
fluid. Boussinesq equations (1 are integrated by a standard 
fully dealiased pseudo-spectral code at resolution N x X N y x N z 
with N y = N x and aspect ratio L x /L z — N x /N z = r (here 
N x = 1024 and r = 1). Other parameters are fig = 0.5, 6o — 1 
(Ag = 0.25), Pr = v/n = 1 and v is chosen such that fc max ?7 > 
1.2 in all runs at final time. Initial perturbation is seeded by 
adding a 10% of white noise to the initial temperature profile 
in a small layer around 2 = 0. 



In the mixing layer turbulent kinetic energy E = 
(1/2) (v 2 ) is produced at the expense of potential energy 
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P = —j3g{zT) as the energy balance indicates 

where e v = v((d a vp) 2 ) is the viscous energy dissipation 
and () represents the integral over the physical domain. 
Assuming that in the turbulent state all quantities in 
© scale in the same way one can balance dv 2 ms /dt ~ 
(3g0QV rms (because temperature fluctuations are bounded 
by the initial jump 9q) and therefore one obtains the 
temporal scaling of velocity fluctuations v rms ~ PgO^t ~ 
Agt, i.e. a motion forced with constant acceleration g. 

The accelerated growth of the width of the mixing 
layer is one of the standard diagnostics in the stud- 
ies of RT turbulence [3, ll4Hl6l |. Several definitions for 
the width have been proposed, based on either local 
or global properties of the mean temperature profile 
T{z 1 t) = l/(L x L y ) J T(x,t)dxdy. The simplest mea- 
sure h r is based on the threshold value of z at which 
T(z,t) reaches a fraction r of the maximum value i.e. 
T(±h r (t)/2,t) = Tr0 /2 Q. This local definition of h 
can be rather noisy and therefore alternative definitions 
based on integral quantities have been proposed 0, HI [ItJ 

c L./2 

h M = I M(c)dz (4) 

J-L z /2 

where c = (T max - T)/(T max - T min ) = 1/2 - Tj 9 is the 
normalized dimensionless temperature (0 < c < 1) and 
M is a mixing function which has support on the mi xing 
layer only, e.g. a logistic function M(c) = 4c(l — c) [10 ] 
or a tent function M(c) = 2c+ (2 - 4c)0(c- 1/2) Di- 
mensionally, h is expected to grow with accelerated law 
h(t) = aAgt 2 with the dimensionless coefficient a which 
depends on the definition of h and apparently also on the 
form of the initial perturbation of the interface 3, ljj . 
Recent studies [!, [3| have shown that a more robust and 
consistent determination of a can be obtained if an initial 
time to is taken into account (physically represent- 
ing the offset at which the t 2 law sets in) suggesting the 
possibility of a universal value, independent on the form 
of the initial perturbation. 

The evolution equation for the normalized tempera- 
ture profile c(z,t) is obtained by averaging ([5} over the 
horizontal directions (assumed periodic) 

d t c + d z wc = nd 2 c (5) 

where w represent the vertical velocity. The thermal flux 
term wc makes (O not closed. Following a common ap- 
proach in turbulence, we close this equation in terms of 
an eddy diffusivity K(z,t) so that (0 is rewritten as 

d t c = d z K(z,t)d z c (6) 

Molecular diffusivity k, included additively in K(z,t), 
can be neglected for large scale properties at high Peclet 
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FIG. 2: (color online). Normalized mean temperature profile 
c(z) computed by averaging over horizontal planes the tur- 
bulent temperature field of Fig. [T] Blue dotted line is the 
prediction of the linear diffusion model 0. Black contin- 
uous line is the fit with the nonlinear model (|10[l. Lower 
inset: enlargement of the temperature profile at the edge of 
the mixing layer. Upper inset: evolution of zi obtained by 
fitting the temperature profile with (I10|l at different times 
and over four different realizations. The line represents the 
fit z\ = jAg(t + to) 2 which gives 7 ~ 0.025 and to — 3.3. 

number. The simplest approximation is to consider K 
independent on z. For our problem, being a diffusion co- 
efficient (i.e. a velocity time a scale) the eddy diffusivity 
is expected to depend on t as K{£) — b 2 (Ag) 2 t 3 with b 
a free dimensionless parameter. The self-similar solution 
to ([5]) with a step initial condition c(z, 0) = 0(z) is 
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FIG. 3: (color online). Heat flux profile wT obtained at the 
same time of Fig. [2] Black line represents the prediction of 
the nonlinear diffusion model as discussed in the text. 

The constant diffusivity solution ([7]) is a relatively good 
approximation of the actual profile obtained from the 
numerical simulations of the full set of equations ([T][2]), 
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as shown in Fig. [2j A closer inspection of the figure 
reveals that the model profile ([7]) is smoother than the 
actual profile at the edges of the mixing region (see inset 
of Fig. The physical origin of this discrepancy is that 
turbulent mixing is not homogeneous within the mixing 
layer. Indeed turbulent velocity fluctuations decrease at 
the ends of the mixing region, and therefore a constant 
K overestimates the diffusivity in these regions. 

An improved model must therefore take into account 
a z-dependence of the diffusivity. Within the general 
framework of mixing length theory by Prandtl jl|, lioj ]. 
the eddy diffusivity can be written as K{z,t) ~ H 2 d z V 
where H represents a length characteristic of mixing and 
V is the typical velocity fluctuation. Because velocity 
is driven by buoyancy at large scale, from equation (fTJ) 
one can estimate that after a time t the typical velocity 
is V oc figTt and taking H cx h(t) one obtains for the 
eddy diffusivity K(z,t) = a(Ag) 3 t 5 d z c where a is again 
a dimensionless constant to be determined empirically. 
We remark that a similar approach, based on gradient 
dependent diffusivity, has been recently used for success- 
fully modeling mixing in stratified flows [2l|. Inserting 
the above expression in one obtains a nonlinear dif- 
fusive model for the mean temperature profile 



d t c = a(Ag) 3 t b d z (d z c) 



(8) 



Observe that the non-linearity of ([8]) reflects the fact that 
temperature fluctuations are not passive in this problem 
as they drive velocity fluctuations in (fTJ). 

Introducing the concentration derivative (p(z,t) = 
3/(a(Ag) 3 )d z c(z,t) and a new time variable t' = t 6 , (j8]) 
is rewritten in a more standard form 



dt><p = d z (ip n d z ip) 



(9) 



with n = 1. Equation ([9]) represents a class of non- 
linear diffusion equations with concentration dependent 
diffusivity well studied in different fields such as thermal 
waves in plasma radiation |22j | and diffusion problems in 
porous media where for our case n = 1 equation @ is 
also known with the name of Boussinesq equation [23|. 
The value of n governs the behavior of the gradient when 
ip — >• which is finite for the present case. The self-similar 
solution (for general n and dimensionality) is known [24j 
and gives for our case 



c(z,t) 



kJL 
4 z. 



c(z,t) = 
c(z,t) = 1 



(*)" 



1*1 < zi 
z < —z\ 

z > Z\ 



(10) 



where z x (t) = ^Agt 2 with 7 = (Sa/2) 1 / 3 . 

Having the analytical expression (fT0|) for the mean con- 
centration, the different definitions of the width of the 
mixing layer are all expressed in terms of z\ and differ 
by a factor only (e.g. h\ = 1z\ and hu = (3/4) zi for 



the tent function Q). Figure [2] shows that the polyno- 
mial function (110[) fits very well the mean concentration 
profile obtained from numerical simulations. Runs at dif- 
ferent resolutions (and viscosity, the only parameter in 
dUt^J when Pr = 1) give analogous results. By fitting 
the numerical profiles at different times, one obtains the 
evolution of z\ displayed in the inset of Fig. [2] which is 
consistent with the quadratic law z\ = yAg(t + t n ) 2 (to 
is the reference time as discussed above). The value ob- 
tained in this way for the coefficient is 7 = 0.025 ± 0.002 
which for the profile Km gives a — (3/4)7 ~ 0.019 in 
agreement with previous numerical results 0, [l(J ■ 

The nonlinear diffusion model can be extended from 
geometrical quantities to study the evolution of dynam- 
ical properties of turbulent convection. In particular, in 
the limit of small thermal diffusivity, from ([5]) and one 
has an expression for the turbulent heat flux in terms of 
the mean temperature profile wT — a(Ag) 3 t 5 (d z T) 2 / 0q. 
Figure [3] shows that the numerically measured profile of 
the heat flux is indeed quite close to the model predic- 
tion, a justification a posteriori of the proposed nonlin- 
ear closure scheme. Using the definition in Q the loss 
of potential energy in kinetic energy (and dissipation) is 
written as —dP/dt = (4/5)j 2 (Ag) 3 t 3 which shows that 
7 is a measure of the efficiency of conversion of available 
potential energy in the turbulent flow. 




FIG. 4: Nusselt number Nu — {wT)/(k9o) versus Rayleigh 
number Ra = Agh 3 / (vk) from three different set of simula- 
tions at resolutions 256 x 256 x 1024 (squares) , 512 x 512 x 2048 
(circles) and 1024 x 1024 x 1024 (triangles) at Pr = 1. 
Kinematic viscosities for the three runs are respectively v = 
6 x 10~ 4 , v = 3 X 10~ 4 and v = 1 x 10" 4 . The line is the 
prediction (fTTf) with 7 = 0.025. 

The relation between the heat flux and the profile ge- 
ometry can be reformulated in terms of dimensionless 
quantities. Indeed, integrating over the width of the mix- 
ing layer it gives a relation between the Nusselt number 
Nu = (wT) / (k9o) (the ratio of convective to conductive 
heat transfer) and the Rayleigh number Ra = Agh 3 /(uk) 
(the ratio of the buoyancy forces to diffusivities). Using 
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the expression (|T0|) and for the length h = hi = 2zi, one 
obtains the temporal evolution laws for the two quanti- 
ties as Ra = 8j 3 {Ag)H e /{vK), Nu = 2 1 2 {Ag) 2 t 3 /(5k) 
and therefore the relation 



Nu 



1 



bV2 



1 ' 2 Pr 1 l 2 Ra>l 2 



(11) 



Equation (|TTj) represents the well known Kraichnan's 
prediction for the "ultimate state of thermal convection" 



25l |26| | which is a regime of turbulent convection ex- 



pected to hold when the contribution of thermal and 
kinetic boundary layers becomes negligible. Because of 
the absence of boundaries, RT turbulence is a natural 
candidate for the appearance of this regime which has 
indeed been observed recently in numerical simulations 
both in two and three dimensions llj, iTJ, |27j . Figure [4] 
shows that the prediction (ITT]) with 7 = 0.025 fits well 
the numerical data obtained from a set of simulations at 
different resolutions. The fact that Nu 3> 1 is a posteri- 
ori confirmation of the negligible contribution of thermal 
diffusivity. 

It is interesting to observe that the above result for 
Nu satisfies a general bound which can be easily ob- 
tained starting from ©• Neglecting thermal diffusiv- 
ity and assuming a self-similar evolution of the pro- 
file c(z,t) = f(z/z\(t)) with the symmetry condition 
f(—z) = 1 — f(z), integrating ([5]) over the z domain 
[— L z /2, L z /2] twice, one obtains 



dzwc 



2ii 



Zl 



dzzc(z, t) 



(12) 



Using the fact that for z > c(z,t) > 1/2 and assuming 
that the flow is still unmixed, c(z, t) = 1 for z > zi, we 
get a bound 



Nu 



1 f— , 1 . 

dzwc < —Z1Z1 

K 



K 



(13) 



If we now further assume the accelerated growth of the 
mixing layer, Zi(t) = jAgt 2 we end with a bound on the 
growth of the Nusselt number 



Nu < -j 2 (Ag) 2 t 3 

K 



(14) 



which is indeed satisfied by our model. The physical in- 
terpretation of this bound is transparent: the growth of 
the heat flux follows the dimensional t 3 law with a coef- 
ficient which depends on the shape of the mean temper- 
ature profile. Maximum growth (|14[) is achieved when 
c(z,t) — 1/2 for —zi < z < zi which means a perfect 
mixing within the mixing layer. This would correspond 
to a coefficient (7/2) 1 / 2 in ([TTj) . 

In this Letter we have introduced a nonlinear diffu- 
sion model with a gradient dependent eddy diffusivity 



which reproduces accurately the large scale phenomenol- 
ogy of Rayleigh- Taylor turbulence obtained from high- 
resolution numerical simulations. The model contains a 
single free parameter, a measure of the turbulence pro- 
duction efficiency, which is directly related to the rate 
of accelerated growth of the mixing layer. The proposed 
closure scheme represents an important step for a phe- 
nomenological description of RT turbulence as it connects 
the evolution of the Nusselt number to the growth of the 
mixing layer, a global geometrical quantity which can be 
easily obtained in experiments. 
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